# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(sdmTMB)
library(devtools)
library(data.table)
library(here)
library(stringr)
library(boot)
library(sdmTMB)
library(patchwork)

# Source code for map plots
# change to url once we have the final one
source("/Users/maxlindmark/Dropbox/Max work/R/marine-litter/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5

# Source utm function
source_url("https://raw.githubusercontent.com/maxlindmark/marine-litter/main/R/functions/lon-lat-utm.R")

# Continuous colors
options(ggplot2.continuous.colour = "viridis")

# Discrete colors
scale_colour_discrete <- function(...) {
  scale_colour_brewer(palette = "Dark2")
}

scale_fill_discrete <- function(...) {
  scale_fill_brewer(palette = "Dark2")
}

Read data

# Read and make data long so that I can for loop through all categories
biom <- read_csv("data/west_coast_litter_biomass.csv") %>% 
  mutate(year_f = as.factor(year),
         quarter_f = as.factor(quarter),
         depth_sc = (depth - mean(depth)) / sd(depth)) %>% 
  rename(plastic = plast) %>% 
  pivot_longer(c(plastic, sup, fishery),
               names_to = "model", # change this to model from litter category for plotting
               values_to = "density") %>% 
  filter(model %in% c("fishery", "plastic", "sup")) %>% 
  mutate(model = str_to_title(model))
#> Rows: 609 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (4): haul_id, ansse_area, trend_area, var
#> dbl (18): diverse, naturmaterial, plast, metall, glas, gummi, quarter, st_no...
#> lgl  (1): balse_area
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: new variable 'year_f' (factor) with 10 unique values and 0% NA
#> 
#>         new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> 
#>         new variable 'depth_sc' (double) with 125 unique values and 0% NA
#> 
#> rename: renamed one variable (plastic)
#> 
#> pivot_longer: reorganized (plastic, sup, fishery) into (model, density) [was 609x26, now 1827x25]
#> 
#> filter: no rows removed
#> 
#> mutate: changed 1,827 values (100%) of 'model' (0 new NA)

num <- read_csv("data/west_coast_litter_numbers.csv") %>% 
  mutate(year_f = as.factor(year),
         quarter_f = as.factor(quarter),
         depth_sc = (depth - mean(depth)) /sd(depth)) %>% 
  rename(plastic = plast) %>% 
  pivot_longer(c(plastic, sup, fishery),
               names_to = "model", # change this to model from litter category for plotting
               values_to = "density") %>% 
  mutate(number = density * swept_area_km2) %>% 
  filter(model %in% c("fishery", "plastic", "sup")) %>% 
  mutate(model = str_to_title(model))
#> Rows: 609 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (4): haul_id, ansse_area, trend_area, var
#> dbl (18): diverse, naturmaterial, plast, metall, glas, gummi, quarter, st_no...
#> lgl  (1): balse_area
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: new variable 'year_f' (factor) with 10 unique values and 0% NA
#> 
#>         new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> 
#>         new variable 'depth_sc' (double) with 125 unique values and 0% NA
#> 
#> rename: renamed one variable (plastic)
#> 
#> pivot_longer: reorganized (plastic, sup, fishery) into (model, density) [was 609x26, now 1827x25]
#> 
#> mutate: new variable 'number' (double) with 125 unique values and 0% NA
#> 
#> filter: no rows removed
#> 
#> mutate: changed 1,827 values (100%) of 'model' (0 new NA)

binom <- read_csv("data/west_coast_litter_biomass.csv") %>% 
  mutate(year_f = as.factor(year),
         quarter_f = as.factor(quarter),
         depth_sc = (depth - mean(depth)) / sd(depth)) %>% 
  rename(other = diverse,
         natural = naturmaterial,
         glass = glas,
         metal = metall,
         rubber = gummi) %>% 
  pivot_longer(c(other, natural, glass, metal, rubber),
               names_to = "model", # change this to model from litter category for plotting
               values_to = "density") %>% 
  filter(!model %in% c("fishery", "plastic", "sup")) %>% 
  mutate(present = ifelse(density == 0, 0, 1)) %>% 
  mutate(model = str_to_title(model))
#> Rows: 609 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (4): haul_id, ansse_area, trend_area, var
#> dbl (18): diverse, naturmaterial, plast, metall, glas, gummi, quarter, st_no...
#> lgl  (1): balse_area
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: new variable 'year_f' (factor) with 10 unique values and 0% NA
#> 
#>         new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> 
#>         new variable 'depth_sc' (double) with 125 unique values and 0% NA
#> 
#> rename: renamed 5 variables (other, natural, metal, glass, rubber)
#> 
#> pivot_longer: reorganized (other, natural, metal, glass, rubber) into (model, density) [was 609x26, now 3045x23]
#> 
#> filter: no rows removed
#> 
#> mutate: new variable 'present' (double) with 2 unique values and 0% NA
#> 
#> mutate: changed 3,045 values (100%) of 'model' (0 new NA)

# Load pred grid
# pred_grid_west <- read_csv("data/pred_grid_west.csv") %>% 
#   mutate(year_f = as.factor(year),
#          quarter_f = as.factor(1),
#          depth_sc = (depth - mean(num$depth)) / sd(num$depth)) %>% 
#   mutate(X = X*1000,
#          Y = Y*1000)

Read model files

# Example

# Read model files
file_paths <- list.files(path = here::here("output/models"),
                         pattern = ".rds",
                         full.names = TRUE)

file_names <- gsub(pattern = ".rds$",
                   replacement = "",
                   x = basename(file_paths))

for(i in 1:length(file_names)){

  assign(file_names[i], readRDS(file_paths[i]))

}

Read csv files with predictions, sims and coefficients

pred_binom <- read_csv("output/dat_pred_binom.csv") %>% mutate(model = str_to_title(model))
#> Rows: 212900 Columns: 11
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (1): model
#> dbl (10): X, Y, year, depth, year_f, quarter_f, depth_sc, est, est_non_rf, e...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 212,900 values (100%) of 'model' (0 new NA)
sim_binom <- read_csv("output/dat_sim_binom.csv") %>% mutate(model = str_to_title(model))
#> Rows: 50 Columns: 7
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (6): year, est, lwr, upr, log_est, se
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 50 values (100%) of 'model' (0 new NA)
sims_binom <- read_csv("output/dat_sims_binom.csv") %>% mutate(model = str_to_title(model))
#> Rows: 25000 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (3): year, .value, .iteration
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 25,000 values (100%) of 'model' (0 new NA)
coef_binom <- read_csv("output/dat_coef_binom.csv") %>% mutate(model = str_to_title(model))
#> Rows: 10 Columns: 6
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): term, model
#> dbl (4): estimate, std.error, conf.low, conf.high
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 10 values (100%) of 'model' (0 new NA)

pred_num <- read_csv("output/dat_pred_num.csv") %>% mutate(model = str_to_title(model))
#> Rows: 127740 Columns: 12
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (1): model
#> dbl (11): X, Y, year, depth, year_f, quarter_f, depth_sc, est, est_non_rf, e...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 127,740 values (100%) of 'model' (0 new NA)
sim_num <- read_csv("output/dat_sim_num.csv") %>% mutate(model = str_to_title(model))
#> Rows: 30 Columns: 7
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (6): year, est, lwr, upr, log_est, se
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 30 values (100%) of 'model' (0 new NA)
sims_num <- read_csv("output/dat_sims_num.csv") %>% mutate(model = str_to_title(model))
#> Rows: 15000 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (3): year, .value, .iteration
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 15,000 values (100%) of 'model' (0 new NA)
coef_num <- read_csv("output/dat_coef_num.csv") %>% mutate(model = str_to_title(model))
#> Rows: 33 Columns: 6
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): term, model
#> dbl (4): estimate, std.error, conf.low, conf.high
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 33 values (100%) of 'model' (0 new NA)

pred_biomass <- read_csv("output/dat_pred_biomass.csv") %>% mutate(model = str_to_title(model))
#> Rows: 127740 Columns: 12
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (1): model
#> dbl (11): X, Y, year, depth, year_f, quarter_f, depth_sc, est, est_non_rf, e...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 127,740 values (100%) of 'model' (0 new NA)
sim_biomass <- read_csv("output/dat_sim_biomass.csv") %>% mutate(model = str_to_title(model))
#> Rows: 30 Columns: 7
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (6): year, est, lwr, upr, log_est, se
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 30 values (100%) of 'model' (0 new NA)
sims_biomass <- read_csv("output/dat_sims_biomass.csv") %>% mutate(model = str_to_title(model))
#> Rows: 15000 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (1): model
#> dbl (3): year, .value, .iteration
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 15,000 values (100%) of 'model' (0 new NA)
coef_biom <- read_csv("output/dat_coef_biomass.csv") %>% mutate(model = str_to_title(model))
#> Rows: 33 Columns: 6
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): term, model
#> dbl (4): estimate, std.error, conf.low, conf.high
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: changed 33 values (100%) of 'model' (0 new NA)
pal <- brewer.pal(n = 8, name = "Paired")

Plot depth coefficients

coefs <- bind_rows(coef_binom, coef_num, coef_biom) %>%
  filter(term == "depth_sc") %>% 
  separate(model, sep = "_", into = c("Model", "Category")) %>% 
  mutate(Category = str_to_title(Category))
#> filter: removed 65 rows (86%), 11 rows remaining
#> mutate: changed 11 values (100%) of 'Category' (0 new NA)
  
ggplot(coefs, aes(y = estimate, x = Category, ymin = conf.low, ymax = conf.high, color = Model)) + 
  geom_point(size = 2, position = position_dodge(width = 0.2)) +
  geom_errorbar(width = 0, position = position_dodge(width = 0.2)) + 
  geom_hline(yintercept = 0, size = 0.75, color = "gray", linetype = 2) + 
  scale_color_manual(values = pal[c(3, 2, 1)]) +
  labs(y = "Estimate") +
  theme_plot() +
  coord_flip() +
  NULL


ggsave("figures/coefs.pdf", units = "cm", width = 10, height = 20)

Plot binomial predictions

dd <- binom %>% 
  group_by(model, year) %>% 
  summarise(mean_present = mean(present)) %>% 
  filter(!model %in% c("Glass", "Rubber"))
#> group_by: 2 grouping variables (model, year)
#> summarise: now 50 rows and 3 columns, one group variable remaining (model)
#> filter (grouped): removed 20 rows (40%), 30 rows remaining

p1 <- sims_binom %>% 
  filter(!model %in% c("Glass", "Rubber")) %>% 
  filter(.iteration < 26) %>% 
  ggplot() +
  facet_wrap(~model, ncol = 3, scales = "free") +
  geom_ribbon(data = sim_binom %>% filter(!model %in% c("Glass", "Rubber")),
                                          aes(year, est, ymin = lwr, ymax = upr), alpha = 0.2) +
  geom_line(aes(year, .value, group = .iteration), alpha = 0.2) +
  geom_line(data = dd, aes(year, mean_present), color = pal[8], size = 0.75, linetype = 2) + # add data
  geom_line(data = sim_binom %>% filter(!model %in% c("Glass", "Rubber")),
            aes(year, est), size = 0.75) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) + 
  theme_plot() +
  theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
  xlab('Year') +
  ylab('Presence prob. in tow') + 
  NULL
#> filter: removed 10,000 rows (40%), 15,000 rows remaining
#> filter: removed 14,250 rows (95%), 750 rows remaining
#> filter: removed 20 rows (40%), 30 rows remaining
#> filter: removed 20 rows (40%), 30 rows remaining

year_diff <- sims_binom %>% 
  filter(!model %in% c("Glass", "Rubber")) %>% 
  filter(year %in% c(min(year), max(year))) %>% 
  pivot_wider(names_from = year, values_from = .value) %>% 
  mutate("2022-2013" = `2022` - `2013`) %>% 
  pivot_longer(c(`2013`, `2022`, `2022-2013`)) %>% 
  ungroup()
#> filter: removed 10,000 rows (40%), 15,000 rows remaining
#> filter: removed 12,000 rows (80%), 3,000 rows remaining
#> pivot_wider: reorganized (year, .value) into (2013, 2022) [was 3000x4, now 1500x4]
#> mutate: new variable '2022-2013' (double) with 1,500 unique values and 0% NA
#> pivot_longer: reorganized (2013, 2022, 2022-2013) into (name, value) [was 1500x5, now 4500x4]
#> ungroup: no grouping variables

diff_sum <- year_diff %>% 
  filter(name == "2022-2013") %>% 
  group_by(model) %>% 
  summarise(diff_lwr = round(quantile(value, probs = 0.025), digits = 3),
            diff_med = round(quantile(value, probs = 0.5), digits = 3),
            diff_upr = round(quantile(value, probs = 0.975), digits = 3)) %>% 
  mutate(label = paste("2022-2013 = ", diff_med, " [", diff_lwr, ",", diff_upr, "]", sep = "")) 
#> filter: removed 3,000 rows (67%), 1,500 rows remaining
#> group_by: one grouping variable (model)
#> summarise: now 3 rows and 4 columns, ungrouped
#> mutate: new variable 'label' (character) with 3 unique values and 0% NA

# Left join so that we can color "sig"
diff_sum <- diff_sum %>% mutate(sig = ifelse(diff_lwr > 0 & diff_upr > 0, "sig.", "not sig."),
                                sig = ifelse(diff_lwr < 0 & diff_upr < 0, "sig.", sig))
#> mutate: new variable 'sig' (character) with one unique value and 0% NA

year_diff <- left_join(year_diff, diff_sum)
#> Joining, by = "model"
#> left_join: added 5 columns (diff_lwr, diff_med, diff_upr, label, sig)
#>            > rows only in x       0
#>            > rows only in y  (    0)
#>            > matched rows     4,500
#>            >                 =======
#>            > rows total       4,500

p2 <- ggplot(data = filter(year_diff, !name == "2022-2013"), aes(name, value)) +
  geom_violin(alpha = 0.5, fill = "grey30", color = NA) +
  geom_jitter(height = 0, width = 0.05, alpha = 0.2, color = "grey30", size = 0.1) +
  geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
               outlier.shape = NA) +
  theme_plot() +
  theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
  facet_wrap(~ model, scales = "free", ncol = 3) +
  xlab('Year') +
  ylab('Presence prob. in tow') + 
  NULL
#> filter: removed 1,500 rows (33%), 3,000 rows remaining

p3 <- ggplot(data = filter(year_diff, name == "2022-2013"), aes(name, value, fill = sig)) +
  geom_violin(alpha = 0.4, color = NA) +
  geom_jitter(height = 0, width = 0.05, alpha = 0.3, color = "grey30", size = 0.2) +
  geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
               outlier.shape = NA) +
  scale_fill_manual(values = pal[c(2, 6)], name = "") +
  theme_plot() +
  theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
  facet_wrap(~ model, scales = "free", ncol = 3) +
  geom_text(data = diff_sum, aes(label = label), inherit.aes = FALSE,
            x = -Inf, y = Inf, size = 2.5, hjust = -0.05, vjust = 2) +
  xlab('') +
  ylab('Diff. in presence prob. per tow') + 
  geom_hline(yintercept = 0, linetype = 2, size = 0.4) +
  NULL
#> filter: removed 3,000 rows (67%), 1,500 rows remaining

p1 / p2 / p3 + plot_annotation(tag_levels = "A")


ggsave("figures/binom_index.pdf", units = "cm", width = 20, height = 20)

Plot Poisson predictions

# Space
pred_num
#> # A tibble: 127,740 × 12
#>         X      Y  year depth year_f quart…¹ depth…²   est est_n…³ est_rf omega_s
#>     <dbl>  <dbl> <dbl> <dbl>  <dbl>   <dbl>   <dbl> <dbl>   <dbl>  <dbl>   <dbl>
#>  1 3.15e5 6.25e6  2013  228    2013       1   0.693 -1.38   -1.38      0       0
#>  2 3.17e5 6.25e6  2013  228    2013       1   0.693 -1.38   -1.38      0       0
#>  3 3.19e5 6.25e6  2013  227    2013       1   0.675 -1.38   -1.38      0       0
#>  4 3.21e5 6.25e6  2013  228.   2013       1   0.702 -1.38   -1.38      0       0
#>  5 3.23e5 6.25e6  2013  226    2013       1   0.656 -1.38   -1.38      0       0
#>  6 3.25e5 6.25e6  2013  229.   2013       1   0.717 -1.38   -1.38      0       0
#>  7 3.27e5 6.25e6  2013  229.   2013       1   0.711 -1.38   -1.38      0       0
#>  8 3.29e5 6.25e6  2013  228.   2013       1   0.686 -1.38   -1.38      0       0
#>  9 3.31e5 6.25e6  2013  229    2013       1   0.711 -1.38   -1.38      0       0
#> 10 3.33e5 6.25e6  2013  229    2013       1   0.711 -1.38   -1.38      0       0
#> # … with 127,730 more rows, 1 more variable: model <chr>, and abbreviated
#> #   variable names ¹​quarter_f, ²​depth_sc, ³​est_non_rf

pred_num2 <- pred_num #%>%
  #filter(year %in% c(2014, 2016, 2018, 2020, 2022))

plot_map_west + 
  geom_raster(data = pred_num2, aes(X, Y, fill = exp(est))) + 
  scale_fill_viridis_c(
    trans = "sqrt",
    name = "Count",
    # trim extreme high values to make spatial variation more visible
    na.value = "yellow", limits = c(0, quantile(exp(pred_num2$est), 0.99))
    ) +
  facet_grid(model~year) +
  ggtitle("Prediction (fixed effects + all random effects)",
          subtitle = paste("Maximum estimated biomass density =", round(max(exp(pred_num2$est)), digits = 3))) +
  geom_sf(size = 0.1)


ggsave("figures/poisson_maps.pdf", units = "cm", width = 20, height = 20)


# Index
dd <- num %>% 
  group_by(model, year) %>% 
  summarise(mean_num = mean(number))
#> group_by: 2 grouping variables (model, year)
#> summarise: now 30 rows and 3 columns, one group variable remaining (model)

p1 <- sims_num %>% 
  filter(.iteration < 26) %>% 
  ggplot() +
  facet_wrap(~model, ncol = 3, scales = "free") +
  geom_ribbon(data = sim_num, aes(year, est, ymin = lwr, ymax = upr), alpha = 0.2) +
  geom_line(aes(year, .value, group = .iteration), alpha = 0.2) +
  geom_line(data = dd, aes(year, mean_num), color = pal[8], size = 0.75, linetype = 2) + # add data
  geom_line(data = sim_num, aes(year, est), size = 0.75) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) + 
  theme_plot() +
  theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
  xlab('Year') +
  ylab('Mean count per tow') +
  NULL
#> filter: removed 14,250 rows (95%), 750 rows remaining

year_diff <- sims_num %>% 
  filter(year %in% c(min(year), max(year))) %>% 
  pivot_wider(names_from = year, values_from = .value) %>% 
  mutate("2022-2013" = `2022` - `2013`) %>% 
  pivot_longer(c(`2013`, `2022`, `2022-2013`)) %>% 
  ungroup()
#> filter: removed 12,000 rows (80%), 3,000 rows remaining
#> pivot_wider: reorganized (year, .value) into (2013, 2022) [was 3000x4, now 1500x4]
#> mutate: new variable '2022-2013' (double) with 1,500 unique values and 0% NA
#> pivot_longer: reorganized (2013, 2022, 2022-2013) into (name, value) [was 1500x5, now 4500x4]
#> ungroup: no grouping variables

diff_sum <- year_diff %>% 
  filter(name == "2022-2013") %>% 
  group_by(model) %>% 
  summarise(diff_lwr = round(quantile(value, probs = 0.025), digits = 2),
            diff_med = round(quantile(value, probs = 0.5), digits = 2),
            diff_upr = round(quantile(value, probs = 0.975), digits = 2)) %>% 
  mutate(label = paste("2022-2013 = ", diff_med, " [", diff_lwr, ",", diff_upr, "]", sep = "")) 
#> filter: removed 3,000 rows (67%), 1,500 rows remaining
#> group_by: one grouping variable (model)
#> summarise: now 3 rows and 4 columns, ungrouped
#> mutate: new variable 'label' (character) with 3 unique values and 0% NA

# Left join so that we can color "sig"
diff_sum <- diff_sum %>% mutate(sig = ifelse(diff_lwr > 0 & diff_upr > 0, "sig.", "not sig."),
                                sig = ifelse(diff_lwr < 0 & diff_upr < 0, "sig.", sig))
#> mutate: new variable 'sig' (character) with 2 unique values and 0% NA

year_diff <- left_join(year_diff, diff_sum)
#> Joining, by = "model"
#> left_join: added 5 columns (diff_lwr, diff_med, diff_upr, label, sig)
#>            > rows only in x       0
#>            > rows only in y  (    0)
#>            > matched rows     4,500
#>            >                 =======
#>            > rows total       4,500

p2 <- ggplot(data = filter(year_diff, !name == "2022-2013"), aes(name, value)) +
  geom_violin(alpha = 0.5, fill = "grey30", color = NA) +
  geom_jitter(height = 0, width = 0.05, alpha = 0.2, color = "grey30", size = 0.1) +
  geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
               outlier.shape = NA) +
  theme_plot() +
  theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
  facet_wrap(~ model, scales = "free", ncol = 3) +
  xlab('Year') +
  ylab('Mean count per tow') + 
  NULL
#> filter: removed 1,500 rows (33%), 3,000 rows remaining

p3 <- ggplot(data = filter(year_diff, name == "2022-2013"), aes(name, value, fill = sig)) +
  geom_violin(alpha = 0.4, color = NA) +
  geom_jitter(height = 0, width = 0.05, alpha = 0.3, color = "grey30", size = 0.2) +
  geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
               outlier.shape = NA) +
  scale_fill_manual(values = pal[c(2, 6)], name = "") +
  theme_plot() +
  theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
  facet_wrap(~ model, scales = "free", ncol = 3) +
  geom_text(data = diff_sum, aes(label = label), inherit.aes = FALSE,
            x = -Inf, y = Inf, size = 2.5, hjust = -0.05, vjust = 2) +
  xlab('') +
  ylab('Diff. in count per tow') + 
  geom_hline(yintercept = 0, linetype = 2, size = 0.4) +
  NULL
#> filter: removed 3,000 rows (67%), 1,500 rows remaining

p1 / p2 / p3 + plot_annotation(tag_levels = "A")


ggsave("figures/poisson_index.pdf", units = "cm", width = 20, height = 20)

Biomass predictions

# Space
pred_biomass
#> # A tibble: 127,740 × 12
#>         X      Y  year depth year_f quart…¹ depth…²   est est_n…³ est_rf omega_s
#>     <dbl>  <dbl> <dbl> <dbl>  <dbl>   <dbl>   <dbl> <dbl>   <dbl>  <dbl>   <dbl>
#>  1 3.15e5 6.25e6  2013  228    2013       1   0.693 -2.60   -2.60      0       0
#>  2 3.17e5 6.25e6  2013  228    2013       1   0.693 -2.60   -2.60      0       0
#>  3 3.19e5 6.25e6  2013  227    2013       1   0.675 -2.60   -2.60      0       0
#>  4 3.21e5 6.25e6  2013  228.   2013       1   0.702 -2.59   -2.59      0       0
#>  5 3.23e5 6.25e6  2013  226    2013       1   0.656 -2.60   -2.60      0       0
#>  6 3.25e5 6.25e6  2013  229.   2013       1   0.717 -2.59   -2.59      0       0
#>  7 3.27e5 6.25e6  2013  229.   2013       1   0.711 -2.59   -2.59      0       0
#>  8 3.29e5 6.25e6  2013  228.   2013       1   0.686 -2.60   -2.60      0       0
#>  9 3.31e5 6.25e6  2013  229    2013       1   0.711 -2.59   -2.59      0       0
#> 10 3.33e5 6.25e6  2013  229    2013       1   0.711 -2.59   -2.59      0       0
#> # … with 127,730 more rows, 1 more variable: model <chr>, and abbreviated
#> #   variable names ¹​quarter_f, ²​depth_sc, ³​est_non_rf

pred_biomass2 <- pred_biomass #%>%
  #filter(year %in% c(2014, 2016, 2018, 2020, 2022))

plot_map_west + 
  geom_raster(data = pred_biomass, aes(X, Y, fill = exp(est))) + 
  scale_fill_viridis_c(
    trans = "sqrt",
    name = "Biomass density [kg]",
    # trim extreme high values to make spatial variation more visible
    na.value = "yellow", limits = c(0, quantile(exp(pred_biomass2$est), 0.99))
    ) +
  facet_grid(model~year) +
  ggtitle("Prediction (fixed effects + all random effects)",
          subtitle = paste("Maximum estimated biomass density =", round(max(exp(pred_biomass$est))))) +
  geom_sf(size = 0.1)


ggsave("figures/tweedie_maps.pdf", units = "cm", width = 20, height = 20)


# Index
p1 <- sims_biomass %>% 
  filter(.iteration < 26) %>% 
  ggplot() +
  facet_wrap(~model, ncol = 3, scales = "free") +
  geom_ribbon(data = sim_biomass, aes(year, est, ymin = lwr, ymax = upr), alpha = 0.2) +
  geom_line(aes(year, .value, group = .iteration), alpha = 0.2) +
  geom_line(data = sim_biomass, aes(year, est), size = 0.75) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) + 
  theme_plot() +
  theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
  xlab('Year') +
  ylab('Biomass [kg]') +
  NULL
#> filter: removed 14,250 rows (95%), 750 rows remaining

p1


year_diff <- sims_biomass %>% 
  filter(year %in% c(min(year), max(year))) %>% 
  pivot_wider(names_from = year, values_from = .value) %>% 
  mutate("2022-2013" = `2022` - `2013`) %>% 
  pivot_longer(c(`2013`, `2022`, `2022-2013`)) %>% 
  ungroup()
#> filter: removed 12,000 rows (80%), 3,000 rows remaining
#> pivot_wider: reorganized (year, .value) into (2013, 2022) [was 3000x4, now 1500x4]
#> mutate: new variable '2022-2013' (double) with 1,500 unique values and 0% NA
#> pivot_longer: reorganized (2013, 2022, 2022-2013) into (name, value) [was 1500x5, now 4500x4]
#> ungroup: no grouping variables

diff_sum <- year_diff %>% 
  filter(name == "2022-2013") %>% 
  group_by(model) %>% 
  summarise(diff_lwr = round(quantile(value, probs = 0.025), digits = 2),
            diff_med = round(quantile(value, probs = 0.5), digits = 2),
            diff_upr = round(quantile(value, probs = 0.975), digits = 2)) %>% 
  mutate(label = paste("2022-2013 = ", diff_med, " [", diff_lwr, ",", diff_upr, "]", sep = "")) 
#> filter: removed 3,000 rows (67%), 1,500 rows remaining
#> group_by: one grouping variable (model)
#> summarise: now 3 rows and 4 columns, ungrouped
#> mutate: new variable 'label' (character) with 3 unique values and 0% NA

# Left join so that we can color "sig"
diff_sum <- diff_sum %>% mutate(sig = ifelse(diff_lwr > 0 & diff_upr > 0, "sig.", "not sig."),
                                sig = ifelse(diff_lwr < 0 & diff_upr < 0, "sig.", sig))
#> mutate: new variable 'sig' (character) with one unique value and 0% NA

year_diff <- left_join(year_diff, diff_sum)
#> Joining, by = "model"
#> left_join: added 5 columns (diff_lwr, diff_med, diff_upr, label, sig)
#>            > rows only in x       0
#>            > rows only in y  (    0)
#>            > matched rows     4,500
#>            >                 =======
#>            > rows total       4,500

p2 <- ggplot(data = filter(year_diff, !name == "2022-2013"), aes(name, value)) +
  geom_violin(alpha = 0.5, fill = "grey30", color = NA) +
  geom_jitter(height = 0, width = 0.05, alpha = 0.2, color = "grey30", size = 0.1) +
  geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
               outlier.shape = NA) +
  theme_plot() +
  theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
  facet_wrap(~ model, scales = "free", ncol = 3) +
  xlab('Year') +
  ylab('Biomass [kg]') + 
  NULL
#> filter: removed 1,500 rows (33%), 3,000 rows remaining

p3 <- ggplot(data = filter(year_diff, name == "2022-2013"), aes(name, value, fill = sig)) +
  geom_violin(alpha = 0.4, color = NA) +
  geom_jitter(height = 0, width = 0.05, alpha = 0.3, color = "grey30", size = 0.2) +
  geom_boxplot(fill = NA, width = 0.1, size = 0.5, color = "grey5",
               outlier.shape = NA) +
  scale_fill_manual(values = pal[c(6, 2)], name = "") +
  theme_plot() +
  theme(plot.margin = unit(c(0, 0.4, 0, 0.4), "cm")) +
  facet_wrap(~ model, scales = "free", ncol = 3) +
  geom_text(data = diff_sum, aes(label = label), inherit.aes = FALSE,
            x = -Inf, y = Inf, size = 2.5, hjust = -0.05, vjust = 2) +
  xlab('') +
  ylab('Diff. biomass') + 
  geom_hline(yintercept = 0, linetype = 2, size = 0.4) +
  NULL
#> filter: removed 3,000 rows (67%), 1,500 rows remaining

p1 / p2 / p3 + plot_annotation(tag_levels = "A")


ggsave("figures/tweedie_index.pdf", units = "cm", width = 20, height = 20)
knitr::knit_exit
#> function (append, fully = TRUE) 
#> {
#>     if (missing(append)) 
#>         append = if (out_format(c("latex", "sweave", "listings"))) 
#>             "\\end{document}"
#>         else if (out_format("html")) 
#>             "</body>\n</html>"
#>         else ""
#>     .knitEnv$terminate = append
#>     .knitEnv$terminate_fully = fully
#>     invisible()
#> }
#> <bytecode: 0x7faec3c332a8>
#> <environment: namespace:knitr>